## CPS RR
## Analysis for "Peacekeepers without helmets: How violence shapes local peacebuilding by civilian peacekeepers"
## Author: Hannah Smidt & Allard Duursma

rm(list=ls())

require(mediation)
require(MASS)
require(foreign)
library(dplyr)
require(sp)

#setwd("SET YOUR OWN WORKING DIRECTORY")

data <- read.dta("./2022_11_07_dataFinal_max_agg_MonEv.dta")

data <- subset(data, !(year==2018 & month>=11) )

data$date <- paste0(1, "_", data$month,"_",data$year)
data$date <- as.Date(data$date, tryFormats = "%d_%m_%Y")
data <- subset(data, !is.na(date))

data$id_str <- paste0(data$admin3RefN, data$admin1Name, data$admin2Name, data$admin3Name)
data$id <- as.numeric(as.factor(data$id_str))

data$un_military_l1_ln = log(data$un_military_l1+0.01)
data$un_military_ln = log(data$un_military+0.01)

data$LocalPeacebuild_AnyAss_ln = log(data$LocalPeacebuild_AnyAssistance+0.01)

####################################################################################
####################################################################################
## MAIN ANALYSIS: Hypothesis 2: Mediation analyses: UN military (lag 1 period), ####
##                  violence (lag 2 periods), and peacebuilding (contemp.)     #####
## Estimations for Figure 8 in Main Manuscript #####################################
####################################################################################
####################################################################################

## Violence in contemporary month

med.fit1a <- glm(un_military_base_l1  ~ ACLED_viol_any_l2 + 
                  roadDensity + distToCapital + Shape_Area +
                  Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + ACLED_viol_any_3m_MA_l2 
                , family = binomial("probit"), data = data)
summary(med.fit1a)

out.fit1a <- glm.nb(LocalPeacebuild_AnyAssistance  ~ ACLED_viol_any_l2 + un_military_base_l1 + 
                     roadDensity + distToCapital + Shape_Area +
                     Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                     ACLED_viol_any_3m_MA_l2, data = data)
summary(out.fit1a)


# Calculate mediation effect (increase boot to 1000 for analyses in paper)
set.seed(1234)
med.out1a <- mediate(model.m = med.fit1a,
                    model.y = out.fit1a,
                    sims = 1000, boot = T,
                    treat = "ACLED_viol_any_l2", mediator = "un_military_base_l1",
                    covariates = NULL,
                    conf.level = 0.95,
                    control.value = 0,
                    treat.value = 1 )

summary(med.out1a)

out1a <- as.data.frame( rbind(  c(median(med.out1a$z.avg.sims), quantile(med.out1a$z.avg.sims, c(0.025, 0.975)) )
                              , c(median(med.out1a$d.avg.sims), quantile(med.out1a$d.avg.sims, c(0.025, 0.975)) ) ) )
 

## Mediation effect of un military numbers (lag 1 period)

med.fit1b <- lm(un_military_l1_ln  ~ ACLED_viol_any_l2 + 
                 roadDensity + distToCapital + Shape_Area +
                 Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                 ACLED_viol_any_3m_MA_l2, data = data)
summary(med.fit1b)

set.seed(3421)
out.fit1b <- glm.nb(LocalPeacebuild_AnyAssistance ~ ACLED_viol_any_l2 + un_military_l1_ln +
                     roadDensity + distToCapital + Shape_Area +
                     Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                     ACLED_viol_any_3m_MA_l2 
                   , data = data)
summary(out.fit1b)

# Calculate mediation effect (increase boot to 1000 for analyses in paper)
set.seed(1234)
med.out1b <- mediate(model.m = med.fit1b,
                    model.y = out.fit1b,
                    sims = 1000, boot = T,
                    treat = "ACLED_viol_any_l2", mediator = "un_military_l1_ln",
                    covariates = NULL,
                    conf.level = 0.95,
                    control.value = 0,
                    treat.value = 1 )
summary(med.out1b)

out1b <-  rbind(as.matrix(out1a)
                , c(median(med.out1b$z.avg.sims), quantile(med.out1b$z.avg.sims, c(0.025, 0.975)) )
                , c(median(med.out1b$d.avg.sims), quantile(med.out1b$d.avg.sims, c(0.025, 0.975)) ) )

out1b <- as.data.frame(out1b)

out1b$id_merge = c(1,2,3,4)

colnames(out1b) <- c("effects", "low", "up", "id_merge")
rownames(out1b) <- c("ADE1", "ACME1", "ADE2", "ACME2")

foreign::write.dta(out1b, "./MediationAnalysis/mediation_military_ci.dta")




#####################################################################
#####################################################################
# Sensitivity analyses for main analayses  (in Appendix B3.2) #######
#####################################################################
#####################################################################


##################################################################################
# Sensitivity analyses: UN Military base (1m lag), violence (2m lag)

# Mediator model
med.fit.for.sens1a <- glm(un_military_base_l1  ~ ACLED_viol_any_l2 + 
                          roadDensity + distToCapital + Shape_Area +
                          Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                          ACLED_viol_any_3m_MA_l2, 
                          family = binomial(link="probit"),
                          data = data)
summary(med.fit.for.sens1a)

# Note that the outcome model is LINEAR and not neg. binomial for this sensitivity test
out.fit.for.sens1a <- lm(LocalPeacebuild_AnyAss_ln ~ ACLED_viol_any_l2 + un_military_base_l1 +
                          roadDensity + distToCapital + Shape_Area +
                          Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                          ACLED_viol_any_3m_MA_l2, 
                          data = data)
summary(out.fit.for.sens1a)

controls = list("roadDensity"=mean(data$roadDensity), "distToCapital" = mean(data$distToCapital)
                ,"Shape_Area" = mean(data$Shape_Area), "Pop_density" = mean(data$Pop_density)
                ,"FoodInsecurity" = mean(data$FoodInsecurity), "ethnicFract" = mean(data$ethnicFract)
                ,"MoslemGroupSize" = mean(data$MoslemGroupSize)
                ,"ACLED_viol_any_3m_MA_l2" = mean(data$ACLED_viol_any_3m_MA_l2))

set.seed(1234)
med.out.for.sens1a <- mediate(model.m = med.fit.for.sens1a,
                             model.y = out.fit.for.sens1a,
                             sims = 100, boot = T,
                             #boot.ci.type="bca",
                             treat = "ACLED_viol_any_l2",
                             mediator = "un_military_base_l1",
                             covariates = controls,
                             conf.level = 0.95,
                             control.value = 0,
                             treat.value = 1 )
summary(med.out.for.sens1a)


# Sensitivity test for direct effect
med.sens.direct1a <- medsens(med.out.for.sens1a
                            , rho.by = as.vector(0.1)
                            , effect.type = "direct"
                            , sims = 100) 
summary(med.sens.direct1a)

png("./Figures/Sensitivity_Test_directEffect_UNBase.png", width = 450, height=450, units="px")
par(mfrow=c(1,1), mar =c(5.1, 4.1, 4.1, 2.1) )
plot(med.sens.direct1a, sens.par = "rho", main = ""
     , ylim = c(-0.1, 0.3)
     , xlim = c(-1, 1), axes =F 
     , smooth.ci=T
     , xlab = "Correlation between errors"
     , ylab = "Average direct effect of violence (2m lag)")
axis(1)
axis(2)
dev.off()


# Sensitivity test for indirect effect
set.seed(1234)
med.sens.indirect1a <- medsens(med.out.for.sens1a
                               , rho.by = as.vector(0.1)
                               , effect.type = "indirect"
                               , sims=100) 
summary(med.sens.indirect1a)
sqrt(0.9)
0.9/0.98

png("./Figures/Sensitivity_Test_Mediator_UNBase.png", width = 450, height=450, units="px")
par(mfrow=c(1,1), mar =c(5.1, 6.1, 4.1, 2.1) )
plot(med.sens.indirect1a, sens.par = "rho", main = ""
     , ylim = c(-0.1, 0.3)
     , xlim = c(-1, 1), axes =F 
     , smooth.ci=T
     , xlab = "Correlation between errors"
     , ylab = "Average mediation effect of\n UN military base (1m lag)" )
axis(1)
axis(2)
dev.off()


##################################################################################################
# Sensitivity analyses: UN Military personnel numbers (1m lag), violence (2m lag)

# Mediator model
med.fit.for.sens1b <- lm(un_military_l1_ln  ~ ACLED_viol_any_l2  ,
                          data = data)
summary(med.fit.for.sens1b)

# Note that the outcome model is LINEAR and not neg. binomial for this sensitivity test
out.fit.for.sens1b <- lm(LocalPeacebuild_AnyAss_ln ~ ACLED_viol_any_l2 + un_military_l1_ln ,
                        data = data)
summary(out.fit.for.sens1b)


set.seed(1234)
med.out.for.sens1b <- mediate(model.m = med.fit.for.sens1b,
                             model.y = out.fit.for.sens1b,
                             sims = 100, boot = T,
                             boot.ci.type="bca",
                             treat = "ACLED_viol_any_l2", mediator = "un_military_l1_ln",
                             covariates = NULL,
                             conf.level = 0.95,
                             control.value = 0,
                             treat.value = 1 )
summary(med.out.for.sens1b)


# Sensitivity for the direct effect
med.sens.direct1b <- medsens(med.out.for.sens1b
                             , rho.by = as.vector(0.1)
                             , effect.type = "direct"
                             , sims = 100) 
summary(med.sens.direct1b)
sqrt(0.6)

png("./Figures/Sensitivity_Test_directEffect_UNMil.png", width = 450, height=450, units="px")
par(mfrow=c(1,1), mar =c(5.1, 4.1, 4.1, 2.1) )
plot(med.sens.direct1b, sens.par = "rho", main = ""
     , ylim = c(-1, 1)
     , xlim = c(-1, 1), axes =F 
     , smooth.ci=F
     , xlab = "Correlation between errors (without covariates)"
     , ylab = "Average direct effect of violence")
axis(1)
axis(2)
dev.off()

# Sensitivity tests for direct effect
set.seed(1234)
med.sens.indirect1b <- medsens(med.out.for.sens1b
                              , rho.by = as.vector(0.1)
                              , effect.type = "indirect"
                              , sims=100) 
summary(med.sens.indirect1b)
sqrt(0.4)

png("./Figures/Sensitivity_Test_Mediator_UNMil.png", width = 450, height=450, units="px")
par(mfrow=c(1,1), mar =c(5.1, 4.1, 4.1, 2.1) )
plot(med.sens.indirect1b, sens.par = "rho", main = ""
     , ylim = c(-1, 1)
     , xlim = c(-1, 1), axes =F 
     , smooth.ci=F
     , xlab = "Correlation between errors (without covariates)"
     , ylab = "Average mediation effect of UN military personnel" )
axis(1)
axis(2)
dev.off()




#######################################################################################
#######################################################################################
## ROBUSTNESS TEST: Hypothesis 2: Mediation analyses: UN military, violence ###########
#                   and peacebuilding ALL CONTEMPORANEOUS (in Appendix D)  ############
# Estimations for Figure D1 in Appendix D #############################################
#######################################################################################
#######################################################################################

## UN military presence

med.fit2a <- glm(un_military_base ~ ACLED_viol_any + 
                   roadDensity + distToCapital + Shape_Area +
                   Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                   ACLED_viol_any_3m_MA
                 , family = binomial("probit"), data = data)
summary(med.fit2a)

out.fit2a <- glm.nb(LocalPeacebuild_AnyAssistance  ~ ACLED_viol_any + un_military_base + 
                      roadDensity + distToCapital + Shape_Area +
                      Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                      ACLED_viol_any_3m_MA, data = data)
summary(out.fit2a)


# Calculate mediation effect (increase boot to 1000 for analyses in paper)
set.seed(1234)
med.out2a <- mediate(model.m = med.fit2a,
                     model.y = out.fit2a,
                     sims = 1000, boot = T,
                     treat = "ACLED_viol_any", mediator = "un_military_base",
                     covariates = NULL,
                     conf.level = 0.95,
                     control.value = 0,
                     treat.value = 1 )

summary(med.out2a)

out2a <- as.data.frame( rbind(  c(median(med.out2a$z.avg.sims), quantile(med.out2a$z.avg.sims, c(0.025, 0.975)) )
                                , c(median(med.out2a$d.avg.sims), quantile(med.out2a$d.avg.sims, c(0.025, 0.975)) ) ) )

out2a <- as.data.frame(out2a)

## UN military numbers 

med.fit2b <- lm(un_military_ln  ~ ACLED_viol_any  + 
                  roadDensity + distToCapital + Shape_Area +
                  Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                  ACLED_viol_any_3m_MA_l1, data = data)
summary(med.fit2b)

set.seed(3421)
out.fit2b <- glm.nb(LocalPeacebuild_AnyAssistance ~ ACLED_viol_any + un_military_ln +
                      roadDensity + distToCapital + Shape_Area +
                      Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                      ACLED_viol_any_3m_MA_l1 
                    , data = data)
summary(out.fit2b)

# Calculate mediation effect (increase boot to 1000 for analyses in paper)
set.seed(1234)
med.out2b <- mediate(model.m = med.fit2b,
                     model.y = out.fit2b,
                     sims = 1000, boot = T,
                     treat = "ACLED_viol_any", mediator = "un_military_ln",
                     covariates = NULL,
                     conf.level = 0.95,
                     control.value = 0,
                     treat.value = 1 )
summary(med.out2b)

rm(out2b)
out2b <- as.data.frame( rbind(out2a,  c(median(med.out2b$z.avg.sims), quantile(med.out2b$z.avg.sims, c(0.025, 0.975)) )
                                , c(median(med.out2b$d.avg.sims), quantile(med.out2b$d.avg.sims, c(0.025, 0.975)) ) ) )


out2b <- as.data.frame(out2b)

out2b$id_merge = c(1,2,3,4)

colnames(out2b) <- c("effects", "low", "up", "id_merge")
rownames(out2b) <- c("ADE1", "ACME1", "ADE2", "ACME2")


foreign::write.dta(out2b, "./MediationAnalysis/mediation_military_ci_App.dta")


#####################################################################
#####################################################################
# Sensitivity analyses for ROBUSTNESS TESTS (in Appendix D) #########
#####################################################################
#####################################################################



#####################################################################
# Sensitivity analyses 
# UN Military Base (contemporary) 

med.fit.for.sens2a <- glm(un_military_base  ~ ACLED_viol_any + 
                           roadDensity + distToCapital + Shape_Area +
                           Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                           ACLED_viol_any_3m_MA, 
                           family = binomial(link="probit"), data = data)
summary(med.fit.for.sens2a)

out.fit.for.sens2a <- lm(LocalPeacebuild_AnyAss_ln ~ ACLED_viol_any + un_military_base +
                          roadDensity + distToCapital + Shape_Area +
                          Pop_density + FoodInsecurity + ethnicFract + MoslemGroupSize + 
                          ACLED_viol_any_3m_MA, 
                          data = data)
summary(out.fit.for.sens2a)

controls = list("roadDensity"=mean(data$roadDensity), "distToCapital" = mean(data$distToCapital)
                ,"Shape_Area" = mean(data$Shape_Area), "Pop_density" = mean(data$Pop_density)
                ,"FoodInsecurity" = mean(data$FoodInsecurity), "ethnicFract" = mean(data$ethnicFract)
                ,"MoslemGroupSize" = mean(data$MoslemGroupSize)
                ,"ACLED_viol_any_3m_MA" = mean(data$ACLED_viol_any_3m_MA))

set.seed(1234)
med.out.for.sens2a <- mediate(model.m = med.fit.for.sens2a,
                             model.y = out.fit.for.sens2a,
                             sims = 100, boot = T,
                             treat = "ACLED_viol_any", mediator = "un_military_base",
                             covariates = controls,
                             conf.level = 0.95,
                             control.value = 0,
                             treat.value = 1 )
summary(med.out.for.sens2a)

# Sensitivity check for direct effect
set.seed(1234)
med.sens.direct2a <- medsens(med.out.for.sens2a
                             , rho.by = c(0.1)
                             , effect.type = "direct"
                             , sims = 100) 
summary(med.sens.direct2a)

# Sensitivity check for indirect effect
set.seed(1234)
med.sens.indirect2a <- medsens(med.out.for.sens2a
                              , rho.by = c(0.1)
                              , effect.type = "indirect"
                              , sims=100) 
summary(med.sens.indirect2a)




#####################################################################
# Sensitivity analyses 
# UN Military personnel numbers (contemporary)


med.fit.for.sens2b <- lm(un_military_ln  ~ ACLED_viol_any , 
                           data = data)
summary(med.fit.for.sens2b)

out.fit.for.sens2b <- lm(LocalPeacebuild_AnyAss_ln ~ ACLED_viol_any + un_military_ln, 
                        data = data)
summary(out.fit.for.sens2b)


set.seed(1234)
med.out.for.sens2b <- mediate(model.m = med.fit.for.sens2b,
                             model.y = out.fit.for.sens2b,
                             sims = 10, boot = T,
                             treat = "ACLED_viol_any", mediator = "un_military_ln",
                             covariates = NULL,
                             conf.level = 0.95,
                             control.value = 0,
                             treat.value = 1 )
summary(med.out.for.sens2b)

# Sensitivity checks for direct effect
med.sens.direct2b <- medsens(med.out.for.sens2b
                             , rho.by = as.vector(0.1)
                            , effect.type = "direct"
                            , sims = 10) 
summary(med.sens.direct2b)

# Sensitivity checks for IN direct effect
set.seed(1234)
med.sens.indirect2b <- medsens(med.out.for.sens2b
                               , rho.by = c(0.1)
                               , effect.type = "indirect"
                               , sims=10) 
summary(med.sens.indirect2b)
sqrt(0.5)
